1 はじめに

RをGISとして使う. ポリゴンオブジェクト上に位置するラインオブジェクトの延長を算出する. それぞれ以下のようにする.

  • ポリゴン:土砂災害警戒区域
  • ライン:道路

土砂災害警戒区域上に位置する道路を特定し,その延長を計算することで防災政策を立案するうえでの一つの判断材料を与える.

2 前準備

使用するライブラリの読み込み,国土数値情報からのダウンロード,openstreetmapからのダウンロードを行う.

2.1 ライブラリの読み込み

library(sf)
library(tidyverse)
library(tmap)
library(mapview)

2.2 国土数値情報から

国土数値情報より石川県のデータに関して以下を読み込む.

コードは以下のとおりである.

# 日本語を扱うのでCP932でのエンコードを忘れない
isikawa <- st_read("N03-170101_17_GML/N03-17_17_170101.shp", options = "ENCODING=CP932", stringsAsFactors=F)
# emergency <- st_read("N10-15_17_GML/N10-15_17.shp", options = "ENCODING=CP932", stringsAsFactors=F) 不使用
sediment <- st_read("A33-17_17_GML/A33-17_17_GML/A33-17_17Polygon.shp", options = "ENCODING=CP932", stringsAsFactors=F)

N03_004が市町村名に相当する. ここを七尾市に指定すればよさそうである.

nanao <- isikawa %>%
  filter(N03_004 == "七尾市")

行政区域を七尾市のみに絞ったものができた. このポリゴン上で土砂災害警戒区域を絞り込む.

#emergency_nanao <- st_intersection(nanao, emergency)
sediment_nanao <- st_intersection(nanao, sediment)

そして投影座標系をEPSG2449とする.

nanao <- st_transform(nanao, 2449)
#emergency_nanao <- st_transform(emergency_nanao, 2449)
sediment_nanao <- st_transform(sediment_nanao, 2449)

mapview()を使ってインタラクティブに確認する.

mapview(nanao) +
  mapview(sediment_nanao)

2.3 openstreetmapから

openstreetmapのapiを使用して七尾市の道路情報を取得する. このapiを使えば指定した矩形範囲内の道路情報を得ることができる. ただし対象とする範囲が広すぎる場合は取得まで時間がかかりすぎることがある. 今回は国道や市道などの種別ごとに取得することで,一度に得るオブジェクトの数を減らすことで対処する. 国道や市道などの属性は下のような名称となる.

  • morotway
  • trunk
  • primary
  • secondary
  • tertialy
  • unclassified
  • residential

矩形範囲は緯度経度での指定となる. 七尾市を囲む範囲は(36.964677, 136.772596,37.198499, 137.070257)のようになる.

openstreetmapのapiについてはpythonでのラッパーがあるのでそれを使用する. Rからpythonを呼び出すにはreticulateパッケージをつかう. pythonを使う下準備は次の通りとなる.

reticulate::use_python(Sys.which(names = "python"))
opass <- reticulate::import(module = "overpass")
api = opass$API() # Rからは$で潜っていく.

たくさん作るので,関数にしておく

geojson_as_sf <- function(way_attr){
  # way_attr は文字列
  bound_box <- "36.964677, 136.772596,37.198499, 137.070257"
  query_opass <- str_c('way["highway"="',way_attr,'"](',bound_box,');(._;>;)')
  attr_api <- api$get(query_opass, verbosity = 'geom')
  way_sf_NodeWay <- st_as_sf(st_read(attr_api))
  way_sf <- way_sf_NodeWay %>%
    filter(highway == way_attr)
  way_sf <- st_transform(way_sf, 2449) # 4326から2449にする
}

道路属性ごとに異なるsfオブジェクトを生成する. 道路属性を文字列として格納しgeojson -> sf する.

# やたらと出力が長くなる. message = F だとhtml出力が止まる.
attr_str <- c("motorway", "trunk", "primary", "secondary", "tertiary", "unclassified", "residential")
for (i in attr_str) {
  now_sf <- str_c(i, "_sf")
  assign(now_sf, geojson_as_sf(i))
}

バウンディングボックスは矩形であり,七尾市の外側の道路も取得している. これを七尾市上のみに絞り込む.

for (i in attr_str) {
  now_sf <- str_c(i, "_sf")
  nanao_sf <- str_c(i, "_nanao") # 道路のsfについて 七尾市だけのものをこれに格納する
  assign(nanao_sf, st_intersection(nanao, get(now_sf))) # nanao_sf <- st_intersection(nanao, atr_sf) とするのをまあ
}

最終的には全道路属性が1つのデータフレームにまとまっているsfオブジェクトを生成する. 現時点ではまだ,道路属性ごとに異なるsfオブジェクトである. これを以下の手順でひとまとめにする.

  1. 道路属性個数分だけあるsfオブジェクトをデータフレームとsfcオブジェクトに分解.
  2. 1.で,道路属性個数分のデータフレームができたので,それをひとまとめにする.
  3. 1.で,道路属性個数分のsfcオブジェクトができたので,それをひとまとめにする.
  4. sfをつくる.2.でできた1つのデータフレームと,3.でできた1つのsfcをくっつける

分解作業,すぐおわる.

for(i in attr_str){
  now_nanao <- str_c(i, "_nanao")
  
  # データフレームだけを now_df に格納
  now_df <- str_c(i, "_df")
  assign(now_df, st_set_geometry(get(now_nanao), NULL))
  
  # sfc だけを now_sfc に格納
  now_sfc <- str_c(i, "_sfc")
  assign(now_sfc, st_geometry(get(now_nanao)))
}

bind_rowにより,道路属性ごとのデータフレームを下に付け足す.

# 格納用のデータフレームを作って,それにバインドしていく
Road_attr_df <- data.frame()
for(i in attr_str){
  now_df <- get(str_c(i, "_df")) # ここでget()すれば1行でデータフレームになる
  Road_attr_df <- bind_rows(Road_attr_df, now_df)
}

sfcも全道路属性をひとまとめにする. リストに格納する.

# 空のリストに放り込む
Road_sfc <- list()

for (i in attr_str) {
  now_sfc <- get(str_c(i, "_sfc")) # ここでget()すれば1行で sfc
  Road_sfc <- c(Road_sfc, i = now_sfc)
  
}

sfオブジェクトをつくる.

Road_sf <- st_sf(Road_attr_df, geometry = Road_sfc)
Road_sf <- st_set_crs(Road_sf, 2449) # 本当は上のやつでいっしょにしたい

道路属性ごとに色分けして図化する.

mapview(nanao) +
  mapview(Road_sf,
          zcol = "highway")

3 土砂災害警戒区域上に限定する

highway属性をすべて有する道路オブジェクトRoad_sfが作られた. 土砂災害警戒区域の領域に含まれる道路延長の算出のため,Road_sfをこの領域内だけに限定したオブジェクトをさらに作成する. 土砂災害警戒区域のsfオブジェクトはsediment_nanaoである.

Road_sedi <- st_intersection(sediment_nanao, Road_sf)
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries

bboxをいきなりでなく,sfcでよい??

p1 <- st_point(c(136.951, 37.001))
p2 <- st_point(c(136.980, 37.025))
interest_sfc <- st_sfc(p1, p2, crs = 4326)
interest_sfc <- st_transform(interest_sfc, 2449)
interest_sfc
## Geometry set for 2 features 
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -19194.69 ymin: 111089.9 xmax: -16608.42 ymax: 113747.7
## epsg (SRID):    2449
## proj4string:    +proj=tmerc +lat_0=36 +lon_0=137.1666666666667 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## POINT (-19194.69 111089.9)
## POINT (-16608.42 113747.7)

特定の範囲のみを拡大して表示する.

tm_shape(nanao, bbox = interest_sfc) +
  tm_polygons() +
tm_shape(sediment_nanao) +
  tm_polygons(col = "red") +
tm_shape(Road_sf) +
  tm_lines() +
tm_shape(Road_sedi) +
  tm_lines(lwd = 5)

4 距離を算出する

ラインオブジェクトの長さの算出はst_length()で実施する.

# 全道路延長の算出
Road_sf <- Road_sf %>%
  mutate(road_length = st_length(Road_sf))
# highway属性ごとに集約して算出
Road_sf_sum <- st_set_geometry(Road_sf, NULL) %>%
  group_by(highway) %>%
  summarise(ALL = sum(road_length))

# 土砂災害警戒区域上の道路延長の算出
Road_sedi <- Road_sedi %>%
  mutate(road_length = st_length(Road_sedi))
# highway属性ごとに集約して算出
Road_sedi_sum <- st_set_geometry(Road_sedi, NULL) %>%
  group_by(highway) %>%
  summarise(sedi = sum(road_length))

土砂災害区域に含まれる道路延長は全道路延長の何パーセントに相当するか.

Road_sedi_sum$sedi / Road_sf_sum$ALL * 100
## Units: [1]
## [1]  4.510658  5.864269  8.554392 10.497554  7.256750  8.416051  9.452831

表形式で整理する.

full_join(Road_sedi_sum, Road_sf_sum) %>%
  mutate(ratio = round(sedi/ALL*100,2)) %>%
  knitr::kable()
## Joining, by = "highway"
highway sedi ALL ratio
motorway 1790.554 [m] 39696.07 [m] 4.51 [1]
primary 3166.712 [m] 54000.11 [m] 5.86 [1]
residential 43262.050 [m] 505729.08 [m] 8.55 [1]
secondary 15146.287 [m] 144283.96 [m] 10.50 [1]
tertiary 10698.641 [m] 147430.19 [m] 7.26 [1]
trunk 4234.333 [m] 50312.58 [m] 8.42 [1]
unclassified 33732.717 [m] 356853.07 [m] 9.45 [1]